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ABSTRACT 

This  thesis  presents  a  Fortran  program  that  numerically  solves  the 
steady-state  matrix  Riccati  equation  of  the  quadratic  cost  optimal 
control  problem.  Each  step  of  the  program  is  presented,  analytically 

and  computationally.  The  check  points  incorporated  in  the  program  and 

the  input  parameters  that  can  be  used  to  assure  a  correct  solution  are 

identified  and  discussed.  Difficulties  encountered  when  verifying  the 
program,  and  the  suggested  solutions,  are  also  presented. 
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I.   INTRODUCTION 

The  Fortran  computer  program  presented  in  this  thesis  provides  a 
relatively  quick  and  reliable  means  for  determining  the  unique, 
symmetric,  steady-state  solution  P  of  the  nonlinear  matrix  Riccati 
equation: 

P  =  0  =  -  PA  -  A'P  -  C'C  +  PBR_1B'P  (1) 

occurring  in  optimal  control  theory.  The  remaining  equation  variables 
are  defined  in  the  following  statement  of  the  quadratic  cost  optimal 
control  problem. 

Given  the  linear,  time-invariant,  completely  controllable  system 
defined  by  the  state  equations: 

x  (t)  =  Ax(t)  +  Bu(t)                        (2) 

x  (0)  =  xQ                                (3) 

where: 

x(t)  is  an  n  x  1  state  vector 

u(t)  is  an  m  x  1  unconstrained  control  vector 

A  is  an  n  x  n  matrix 

B  is  an  n  x  m  matrix, 

determine  the  control  vector  u*(«)  which  minimizes  the  quadratic  cost 

functional : 

J(xQ;u(-))  -   /  [  x'(t)C'Cx(t)  +  u'(t)Ru(t)  ]  dt        (4) 
o 

where: 

C  is  an  m  x  n  matrix 

R  is  an  m  x  m  symmetric,  positive  definite  matrix 

and  the  matrix  pair  (A,C)  is  completely  observable. 


It  has  been  shown  [Ref.  1]  that  u*(-)  is  given  by  the  linear  feed- 
back law: 

u*(t)  =  -R"1B'Px(t)  =  -L*x(t)  (5) 

where  P  is  the  unique  solution  to  matrix  equation  (1). 

For  the  system  of  equation  (2)  to  be  completely  controllable,  the 
n  x  mn  augmented  matrix  G, 

G  =  (B  |  AB  |  A2B  |  •••  |  An"]B  ),  (6) 

must  contain  n  linearly  independent  column  vectors,  or  equivalently, 
the  rank  of  G  must  be  n.  For  the  matrix  pair  (A,C)  to  be  completely 
observable,  the  n  x  mn  augmented  matrix  H, 

H  =  (  C  |  A'C  |  (A')2C  |  •••  |  (A')n"1C  ), 

must  contain  n  linearly  independent  column  vectors  or  the  rank  of  H 

must  be  n. 

When  u*(-)  is  given  by  equation  (5),  equation  (2)  can  be  rewritten 

as: 

x(t)  =  Ax(t)  -  BL*x(t)  =  (A  -  BL*)x(t)  (7) 

which,  using  equation  (3),  has  the  general  solution: 

x(t)  -  e<A  "  **>»  x0  .  (8) 

Because  the  system  is  completely  controllable,  as  time  approaches 
infinity  x(t)  must  remain  bounded.  To  satisfy  this  requirement,  the 
matrix  (A  -  Bl_*)  must  be  a  stable  matrix  or,  alternatively,  all  the 
eigenvalves  of  the  matrix  (A  -  BL*)  must  have  negative  real  parts.  As 
a  consequence  of  this: 

!im  x(t)  =  *  lim  e(A  "  BL*>  -  0  .  (9) 


To  evaluate  the  quadratic  cost  associated  with  the  optimal  con- 
trol u*(t),  substitute  equation  (5)  into  equation  (4)  and  write: 

/oo 
[x'(t)C'Cx(t)  +  x'(t)L*'RL*x(t)]  dt.    (10) 

Substituting  equation  (8)  for  x(t),  expanding  L*'RL*,  and  factoring  out 
the  common  terms: 


j*  =  x0 


/-e(A-BL*)'t[c,c  +  pBR-lB,p]e(A-BL*)tdt 
oJ 


v  <"> 


-1D. 


From  equation  (1)  C'C  =  -PA  -A'P  +  PBR  B'P  or: 

C'C  +  PBR'Vp  =  -P(A  -  BL*)  -  (A  -  BL*)'P 


(12) 


Substituting  equation  (12)  into  equation  (11)  and  integrating  by  parts 
the  first  term  of  the  resulting  integral  expression,  we  get: 

re(A-BL*)'tr_p(A_BL*)]e(A-BL*)tdt 
oJ 

(A  -  BL*)'t  [p]  e(A  -  BL*)t 


=  -  e 


t=0 


f   (A-BL*)'  e(A  "  BL*},t  [P]  e 
oJ 


(A  -  BL*)t 


dt  . 


(13) 


Using  equation  (9)  to  evaluate  the  upper  limit  of  the  first  term,  the 
first  term  reduces  to  P.  Recognizing: 


ZeZt  =  Z  I    Jf-zV  =eZtZ 
i=0  1" 


(14) 


we  see  that  the  second  term  in  equation  (13)  will  cancel  the  second 
term  in  equation  (11)  when  equation  (12)  is  substituted.  Thus: 

J*  =  x'Px^  .  (15) 
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II.  ANALYTICAL  METHODS  OF  SOLUTION 

A.  ANALYTICAL  METHOD  OF  KLEINMAN 

The  main  computer  program  of  this  thesis  is  based  upon  a  method  for 
solving  the  steady-state  Riccati  equation  published  by  David  L.  Kleinman 
[Ref.  2]  in  1968.  The  iteration  scheme  for  solving  euqation  (1)  is 
presented  in  the  following  theorem  by  Kleinman. 

1 .  Kleinman 's  Theorem 

Let  V. ,  k  =  0,  1 ,  2,  . . .  be  the  n  x  n  (unique)  positive  definite 
matrix  solution  of  the  linear  algebraic  matrix  equation: 

0  =  A£Vk  +  VkAk  +  C'C  +  L^RLk  (16) 

where,  recursively, 

Lk  =  R"1B,Vk_1        k  =  1,2,3,... 

Ak  -  A  -  BLk         k  =  0,1,2,... 

and  where  L„  is  chosen  such  that  the  matrix  A  =  A  -  BL„  has  eigenvalues 
o  o       o      3 

with  negative  real  parts. 

Then:   1)  P  <  Vk+1  <  Vk  <  • • ■    k  =  0,1,2,... 

2)   lim  v  =  p 

The  notation  X  >  Y  [X  >  Y]  means  that  the  matrix  X  -  Y  is 
positive  [semi]  definite. 

2.  The  Cost  Matrix 

The  proof  of  this  theorem  requires  the  introduction  and  defini- 
tion of  a  cost  matrix  V,  .  Assume  that  u.(x(t))  =  -  Lx(t)  is  an  arbitrary 
feedback  law,  with  feedback  gains  of  L,  and  u. (x(t))  is  applied  to  the 


system  of  equation  (2).  Following  a  development  similar  to  equations 
(7),  (8)  and  (15),  the  resulting  quadratic  cost  function  is: 


J(x0;uL(x(t)))  =  x^VLxo 


where  V,  is  the  cost  matrix  associated  with  the  feedback  gains  L  and  is 
given  by: 

V,  =  f   e(A"BL),t  [C'C  +  L'RL]  e(A  ~  BL)t  dt  .        (17) 

V.  is  bounded  if  and  only  if  the  closed-loop  system  control  matrix  (A  -  BL) 
is  stable.  If  V.  is  bounded,  it  becomes  the  unique  (positive  definite) 
solution  of  the  linear  matrix  equation: 

0  =  (A  -  BL)'VL  +  VL(A  -  BL)  +  C'C  +  L'RL  .  (18) 

Examining  the  first  term  of  equation  (18)  and  substituting 
equation  (17)  for  V.  ,  we  can  verify  this  relationship: 

(A  -  BL)'V  =   /  (A  -  BL)'  e(A  "  BL) ^[C'C  +  L'RL]  e(A  "  BL)t  dt. 

Integrating  by  parts  and  using  equation  (14)  we  have: 

oo 

(A  -  BL)'V1  ■  e(A  ■  BL),t  [C'C  +  L'RL]  e(A  "  BL)t  t=Q 

-  Te(A-  BL>Vc  +  L'RL]e(A-BL)t  (A  -  BL) 
oJ 

Using  equation  (9)  to  evaluate  the  upper  limit  of  the  first  term,  the 

first  term  reduces  to  C'C  +  L'RL  and  we  have, 

(A  -  BL)'VL=-C'C  -  L'RL  -  VL(A  -  BL)  , 

the  desired  result. 


dt 
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Recalling  the  result  of  equation  (15)  we  see  that  for  the 
optimal  control  of  equation  (5)  we  have: 

VL*  =  P  •  (19) 

If  Li  and  L_2  are  the  gain  matrices  associated  with  the  cost 
matrices  V\  and  V2,  it  can  be  shown  [Ref.  3]  that: 

Vl  "  V2  =  J     e(A  "  BL2)t  [(L1  '  L2)'R(L1  "  L2) 
o 

-  (L1  -  L2)'(B'V1  -  RL2) 

-  (B'V1  -  RL2)'(L1  -  L2)]  e(A  "  BL2)t:  dt         (20) 


or: 


V 


1  -  V2  =  J    e(A  "  Bliyt   [(L1  -  L2)'R(L1  -  L>) 


-  (L1  -  L2)'(B'V2  -  RL2) 

-  (B'V2  -  RL2)'(L1  -  L2)]  e(A  '  BLl}t  dt.         (21) 

If  either  matrix  (A  -  Bl_2)  or  matrix  (A  -  BL-. )  is  unstable,  V, 
or  Vp,  respectively,  will  be  unbounded  and  care  must  be  exercised  in 
using  the  above  formulas. 

3.  Proof  of  Kleinman's  Theorem 

1)  P  <  Vk+1  <  Vk  s  ••■    k  =  0,1,2, •••.   Let  VQ  be  the  cost 
matrix  for  L  ,  the  initial  guess  that  yields  a  stable  closed-loop 
system  control  matrix,  and  let  V-,  be  the  cost  matrix  for  L-.  =  R~  B'V  . 
Substituting  V  and  V-,  into  equation  (20)  and  noting: 

eA'kt  =  (eAkt)' 

R  =  Y'Y,  where  Y  is  unique 

B'V0  -  RL-|  =  0 
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we  have: 

Vo  "  Vl  =   /  ^  [(Lo  ■  Ll),R(Lo  "  Ll)]  eA]t  dt  * 

Define  Z(t)  =  Y(LQ  -  L-))eAlt.  It  has  been  shown  [Ref.  4]  that  for 
Z(t)  a  real  matrix: 

Z'(t)Z(t)  ;>  0     for  all  t  *  0. 

/oo 
Z'(t)Z(t)  dt  >  0  or: 

V0  *  V1   .  (22) 

Let  V^  be  the  cost  matrix  associated  with  L*,  use  equation  (19) 
and  substitute  into  equation  (21).  Noting  that: 

B'P  -  RL*  =  0 

we  have: 

V]  -  P  =   /  e^  [(L]  -  L*)'R(L1  -  L*)]  eAlt  dt. 
o 

This  time  define  Z(t)  =  Y(L1  -  L*)eAlt  and  we  have: 


■I 


V1  -  P  =   /  Z'(t)Z(t)  dt  >  0  or: 

V1  ;>  P  .  (23) 

Combining  the  results  of  equations  (22)  and  (23)  we  get: 

meaning  that  V-.  is  bounded  above  by  P  and  below  by  V  .  Thus,  A-|  is  a 
stable  matrix  and  V,  satisfies  equation  (16)  with  k  =  1.  Similar 
arguments  can  be  made  for  k  =  2,3,4,...  yielding: 

the  desired  result. 
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2)  Jlm  Vw  =  P.  The  J1m  V.  =  V  exists  by  a  theorem  on  monotonic 
convergence  of  positive  operators  [Ref.  5].  Thus,  in  the  limit  V|<  = 
Vk+1  and  Ak  =  A  -  BLk  =  A  -  BR~Vvk  =  A  -  BR'Vv^.  Since  A£  = 
A'  -  V|BR~  B'  equation  (16)  becomes,  in  the  limit  as  k  approaches 
infinity, 

0  =  A'Vk  -  VkBR"Vvk  +  VkA  -  VkBR-1B'Vk  +  C'C  +  VkBR""Vvk 

which  is  equation  (1),  the  desired  result. 


B.  BASS'S  METHOD  OF  DETERMINING  A  STABLE  INITIAL  GUESS 

Kleinman's  method  requires  an  initial  guess  of  the  feedback  gain 
that  yields  a  stable  closed-loop  system  control  matrix.  Kleinman 
remarks  [Ref.  2]  that,  if  the  system  of  equation  (2)  is  completely 
controllable,  then  the  desired  initial  guess  will  always  exist.  A 
method  that  could  be  programmed  for  a  general,  controllable  system  to 
yield  this  correct  initial  guess  was  sought. 

Bass  [Ref.  6]  presented,  but  did  not  publish,  a  general  method  for 
determining  a  stable  initial  guess  for  a  completely  controllable  system 
in  1961.  The  method  was  published  in  a  paper  by  Wonham  and  Cashman 
[Ref.  7]  and  a  proof  can  be  found  in  a  subsequent  paper  by  Bass  [Ref. 8]. 
Given  the  controllabe  matrix  pair  (A,B),  the  matrix: 

k     =  A  -  BL„ 
o       o 

will  be  stable  when  L  =  B'X"  ,  where  X  is  the  (unique)  positive  definite 
solution  to  the  linear  matrix  equation: 

-  (A  +  BI)  X  -  X(A  +  61)'  +  2BB'  =  0  (24) 

where  6  is  defined  by: 

3  >  I IAI I 
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where  |-|  is  the  Euclidean  norm.  According  to  Wonham  [Ref.  7]  the 
results  are  also  valid  if: 

, MAX  r  , 

3  >  /IT  .  I        a.. 


J  fa    '  U 


where  a. .  are  the  elements  of  A, 


C.   SOLVING  THE  MATRIX  EQUATION  Y'X  +  XY  =  Z 

Equations  (16)  and  (24)  are  in  the  form  of  the  Lyapunov  equation: 

Y'X  +  XY  =  Z  (25) 

The  methods  found  in  the  literature  for  solving  this  equation  fall  into 
two  general  categories:  series  solutions  and  simultaneous  linear 
equation  solutions. 

In  the  series  solutions,  X  is  found  from  the  sum  of  a  converging 
matrix  series,  i.e.,  Ref.  9.  For  the  series  to  converge  the  matrix  Y 
must  be  a  stable  matrix.   In  solving  equation  (16)  this  condition  is 
met;  A.  is  stable  by  the  definition  of  a  bounded  cost  matrix.   In 
general,  however,  the  matrix  (A  +  pi)'  of  equation  (24)  is  not  stable, 
and  the  series  method  fails  to  solve  equation  (24)  properly. 

Since  the  unique  solution,  X,  is  symmetric,  equation  (25)  represents 
n(n  +  l)/2  unknowns.  The  second  category  of  solutions  expands  equation 
(25)  into  a  set  of  n(n  +  l)/2  simultaneous  linear  equations.  An 
economical  way  of  recursively  expanding  equation  (25)  was  given  by 
Bingulac  [Ref.  10],  using  an  integer  coefficient  matrix  to  expand  the 
equation.  This  is  the  method  used  to  expand  equations  (16)  and  (24) 
to  the  form  Ax  =  B,  which  can  be  solved  by  a  variety  of  simultaneous 
equation  solvers. 
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D.  ALTERNATE  ANALYTICAL  METHODS 

There  are  several  alternate  methods  found  in  the  literature  for 
solving  the  steady-state  matrix  Riccati  equation  (1). 

A  method  developed  by  Bass  [Ref.  11]  obtains  the  solution  from  a 
2n-dimensional  Hamiltonian,  H,  and  the  terms  of  the  polynomial  expansion 
of  the  stable  roots  of  H.  This  scheme  is  considered  too  sensitive  to 
finite  numerical  computations  to  be  of  practical  use. 

MacFarlane  [Ref.  12]  shows  that  the  solution  can  be  obtained  from 
the  eigenvectors  corresponding  to  the  unstable  eigenvalues  of  a 
similar  Hamiltonian.  The  scheme  requires  that  the  system  have  distinct 
eigenvalues  and  that  the  corresponding  eigenvectors  be  determined  (this 
is  difficult  for  high  order  systems). 

Blackburn  [Ref.  13]  programmed  a  method  based  on  a  Newton-Raphson 
iterative  solution  for  simultaneous  nonlinear  equations.  This  scheme 
requires  an  initial  guess  that  yields  a  stable  closed-loop  system  con- 
trol matrix.  The  user  must  determine  this  initial  guess  so  that  it  is 
close  enough  to  the  final  solution  for  the  local  Newton-Raphson  method 
to  converge. 

A  fourth  method  integrates  the  full,  nonsteady-state  Riccati 
equation  (1)  backwards  in  time,  from  a  set  of  zero  initial  conditions, 
until  each  element  of  the  P  matrix  reaches  a  satisfactory,  small  value. 
For  systems  of  even  moderate  order,  this  method  is  prohibitive  with 
respect  to  computation  time. 
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III.  COMPUTATIONAL  METHOD  OF  SOLUTION 

A.  SUBROUTINE  RICATS 

1 .  General 

Subroutine  RICATS  is  the  subroutine  that  the  user  will  call  to 
solve  equation  (1).  The  program  language  is  FORTRAN  IV  for  the  Operat- 
ing System  /  360  which  is  compatible  with,  and  encompasses  USASA 
FORTRAN.  The  calling  arguments,  in  the  required  sequence,  are  explained 
in  the  comment  cards  at  the  beginning  of  the  subroutine.  It  should  be 
noted  that  IA  =  n  and  JB  =  m.  Basically  the  subroutine  iterates  on 
equation  (16),  using  equation  (24)  to  calculate  the  initial  guess,  until 
the  solution  converges. 

For  first  order  systems  (n  =  1),  the  nonlinear  equation  (1) 
becomes  a  quadratic  equation.  For  these  systems,  subroutine  RICATS 
solves  the  resulting  quadratic  and  returns  the  largest  root  in  the 
first  element  of  the  P  matrix  . 

2.  The  System  Controllability  Check 

An  analytical  requirement  of  Bass's  method  of  determining  the 
initial  guess  is  that  the  matrix  pair  (A,B)  be  completely  controllable. 
To  check  this  requirement  the  augmented  matrix  G  of  equation  (6)  is 
formed  and  subroutine  GMRANK  is  called  to  determine  the  rank  of  G.   If 
the  rank  equals  n,  execution  continues,  otherwise  subroutine  RICATS 
returns  with  IER  =  2. 

3.  The  Initial  Guess  Stability  Check 

Kleinman's  iteration  scheme  requires  that  the  eigenvalues  of  the 
closed-loop  system  control  matrix  of  the  initial  guess  have  negative 
real  parts  in  order  for  the  initial  cost  matrix  V  to  be  bounded. 
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In  subroutine  RICATS,  when  the  stability  check  is  requested,  the  matrix 
A  -  BL  is  copied  into  a  work  matrix  and  then  passed  to  the  IBM/SSP 
subroutines  HSBG  and  ATEIG.  The  eigenvalues  computed  by  ATEIG  are  then 
individually  tested  to  be  algebraically  less  than  minus  the  absolute 
value  of  EIGMAX,  where  EIGMAX  is  set  by  the  user.  If  any  eigenvalue 
fails  the  check  RICATS  returns  with  IER  =  3. 

4.  The  Solution  Positive  Definite  Check 

Kleinman's  iteration  scheme  is  supposed  to  converge  to  a  positive 
definite  solution.  The  iterations  can  converge  to  a  nonpositive  definite 
solution  if  all  of  Kleinman's  theorem  requirements  are  not  numerically 
met.  Therefore,  for  the  user's  convenience,  a  check  to  verify  that  the 
solution  is  positive  definite  is  included.  In  subroutine  RICATS,  when 
the  positive  definite  check  is  requested,  the  solution  is  factored  by 
the  Cholesky  square-root  method.  A  nonpositive  term  on  the  diagonal 
of  the  factorized  matrix  constitutes  a  nonpositive  definite  solution 
and  RICATS  returns  with  IER  =  4  +  KL,  where  KL  =  0  is  for  a  converged 
solution  and  KL  =  1  is  for  NTRY  iterations  without  convergence. 

5.  The  Steps  of  Subroutine  RICATS 

Calling  subroutine  RICATS  causes  the  following  sequential 
steps  to  be  executed. 

1.  Check  input  parameters  IA,  JB,  IER,  NTRY  to  insure  they  are 
within  proper  bounds.  If  check  fails  IER  =  6,  NN  =  0  and  RICATS  returns. 

2.  If  the  user  requests,  check  the  controllability  of  the 
inputed  system.  If  check  fails  IER  =  2,  NN  =  0  and  RICATS  returns. 

3.  Set  NN  =  0. 

4.  If  this  is  a  first  order  system  (n  =  1),  go  to  step  30. 

5.  Form  E  =  -  BR"]B'. 
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6.  Form  F  =  2BB\ 

7.  Form  P  =  (A  +  el)'  using  equation  (24)  to  define  3  where 
the  "greater  than"  magnitude  is  provided  by  the  variable  FIX  or: 

3=  FIX^T^X  I     I  a..  |  • 

8.  Solve  (A  +  3l)X  +  X(A  +  3l)'  =  2BB\ 

9.  If  subroutine  SIMQ,  through  subroutine  MLIAPS,  returns  a 
nonzero  error  flag,  NN  =  1. 

10.  Form  L  =  BX"1 . 

o 

11.  Form  P  =  A  -  BLQ  . 

12.  If  the  user  requests,  check  the  stability  of  the  system 
matrix  of  the  initial  guess.  If  the  check  fails,  IER  =  3  and  RICATS 
returns. 

13.  Form  F  =  -  Q  -  L'RL  . 

^    oo 

14.  Solve  (A  -  BLo)'V1  +  V1 (A  -  BLQ)  =  -  Q  -  L^RLQ. 

15.  If  subroutine  SIMQ,  through  subroutine  MLIAPS,  returns  a 
nonzero  error  flag,  NN  =  NN  +  1 . 

16.  Set  k  =  1,  KL  =  0. 

17.  Form  P  =  EVk< 

18.  Form  F  =  -  Q  +  VkP. 

19.  Form  P  =  A  +  P. 

20.  Solve  (A  +  EVk)'Vk+1  +  Vk+1(A  ♦  EVk)  -  -  Q  +  VkEVk. 

21.  If  subroutine  SIMQ,  through  subroutine  MLIAPS,  returns  a 
nonzero  error  flag,  NN  =  NN  +  1 .  If  NN  >  (n  +  l)/2,  go  to  step  29. 

22.  Check  each  element  of  V.+,  by  ABS  (  (v^  -  v.+-i)/vk  ) 
<  TOLER.   If  all  elements  of  V,+,  pass  this  test,  go  to  step  27. 

23.  Form  VR  =  Vk+,  . 
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24.  If  k  >  NTRY;  go  to  step  26. 

25.  Set  k  =  k+1;  go  to  step  17. 

26.  Set  KL  =  1. 

27.  If  the  user  requests,  check  to  see  if  V.+,  is  a  positive 
definite  matrix.  If  check  fails  IER  =  4  +  KL  and  RICATS  returns. 

28.  Set  IER  =  KL  and  RICATS  returns. 

29.  Set  IER  =  7  and  RICATS  returns. 

30.  If  the  user  requests  the  controllability  check  and  B(l)  is 
zero,  IER  =  2,  NN  =  0  and  RICATS  returns. 

31.  Set  P  =  AR/BB  +  ^AR/BB)2  +  C'CR/BB 

32.  If  the  user  requests  the  positive  definite  check  and  P  <, 
0.1E-35,  IER  =  4,  NN  =  0  and  RICATS  returns. 

33.  Set  IER  =  0  and  RICATS  returns. 

B.  SUBROUTINE  MLIAPS 

Subroutine  MLIAPS  is  an  auxiliary  routine  used  by  subroutine  RICATS. 
The  subroutine  expands  the  equation  in  steps  8,  14  and  20  of  subroutine 
RICATS  into  n(n  +  l)/2  simultaneous  linear  equations  of  the  form  Ax  =  B. 
The  method  used  was  presented  by  Bingulac  [Ref.  10]  in  1970.  These 
n(n  =  l)/2  simultaneous  linear  equations  are  then  solved  by  subroutine 
SIMQ.  Upon  return  from  subroutine  SIMQ,  MLIAPS  immediately  returns  to 
subroutine  RICATS  and  any  error  codes  returned  by  SIMQ  are  passed, 
unchanged,  to  RICATS. 

C.  SUBROUTINE  GMRANK 

Subroutine  GMRANK  is  a  general  subroutine  for  determining  the  mini- 
mum row  or  column  rank  of  an  arbitary  real  matrix.  The  method  used  is 
simple  row  and  column  interchanges  for  maximum  pivoting,  with  successive 
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reduction  on  the  remaining  matrix  elements  [Ref.  14].  When  the  absolute 
value  of  a  pivot  element  is  less  than  0.1E-35,  or  when  the  final  pivot 
element  has  been  determined,  the  subroutine  returns  the  integer  K, 
where  K  is  the  number  of  successfully  determined  nonzero  pivot  elements 
or,  equivalently,  the  rank  of  the  inputed  matrix. 
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IV.  OPERATIONAL  ASPECTS  OF  USING  SUBROUTINE  RICATS 

A.  CONVERGENCE  PROPERTIES 

Kleinman  suggests  that  equation  (16)  will  converge  to  a  satisfactory 
solution  within  seven  to  ten  iterations,  using  an  initial  guess  generated 
by  hand  or  from  a  previous  solution.  Subroutine  RICATS,  using  a  correct 
set  of  input  parameters  and  Bass's  method  of  generating  an  initial  guess, 
converged  within  forty  iterations  for  e^ery   system  tested. 

To  test  subroutine  RICATS,  over  seventy-five  systems  were  run  and  a 
satisfactory  solution  was  returned  for  each  one.  As  a  test  of  how 
accurate  the  solutions  were,  the  average  absolute  value  of  the  matrix  P 
of  equation  (1)  was  determined.  For  the  satisfactory  solutions  the 
average  values  were  of  the  order  of  10"  .  The  standard  input  parameters 

were: 

NTRY  =  50 

TOLER  =  0.1E-03 

FIX  =  1.1 

EIGMAX  =  0.001 

IER  =  0 

and  the  usual  error  flag  returned  was  IER  =  0  (see  the  discussion  of 

the  parameter  IER  for  high  order  systems). 

B.  INPUT  PARAMETERS 

Through  its  input  parameters  subroutine  RICATS  was  written  to  be  as 
flexible  as  possible  for  the  user.  For  any  system,  the  five  parameters 
NTRY,  TOLER,  FIX,  EIGMAX,  IER  influence  the  final  returned  solution. 
It  is  hoped  that  the  user  can  tailor  these  input  parameters  to  meet 
his  particular  needs. 

In  the  following  discussion,  "high"  or  "higher  order  systems," 
refer  to  systems  of  the  eighth  order  and  above  (n  >  8). 
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1.  NTRY 

NTRY  is  the  maximum  number  of  iterations  the  user  desires  to 
be  attempted,  before  returning  to  the  user's  program  without  a  converged 
solution.  A  recommended  value  is  fifty  and  the  user  is  usually  assured 
that  a  solution  that  is  going  to  converge,  will  converge  within  fifty 
iterations.  It  should  be  noted  that  some  initial  guesses  generated  may 
yield  marginally  stable  system  control  matrices,  and  these  solutions  will 
require  a  larger  number  of  iterations  to  converge.  From  experience,  one 
hundred  and  fifty  iterations  was  considered  to  be  a  sufficient  practical 
upper  limit. 

2.  TOLER 

TOLER  is  the  maximum  percentage  difference,  between  each  element 
of  the  solution,  on  successive  iterations  for  the  convergence  criterion 
to  be  satisfied.  The  accuracy  of  single-precision  computations  is  about 
five  significant  digits,  so  decreasing  TOLER  beyond  0.0001  will  not 
result  in  an  appreciably  more  accurate  solution.  Decreasing  the  magni- 
tude of  TOLER  (up  to  this  limit)  will  tend  to  increase  the  accuracy 
and  the  number  of  iterations  required  for  the  desired  solution. 

3.  FIX. 

FIX  is  the  constant  used  to  insure  the  magnitude  of  $  in  equa- 
tion (24)  is  greater  than  the  Euclidian  norm  of  the  A  matrix. 
Analytically  then,  FIX  should  be  greater  than  one  and  the  suggested 
value  of  1.1  worked  well  for  the  majority  of  systems.  However,  for 
some  systems  equation  (24)  generates  a  singular  initial  guess,  evidenced 
by  an  underflow  error  message  from  subroutine  MINV.  These  singular 
initial  guesses  can  still  yield  quite  satisfactory  solutions  within 
five  or  six  iterations.  However,  the  user  can,  by  varying  FIX  through 
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the  suggested  range,  remove  this  underflow  from  his  execution.  From 
practical  experience,  the  range  of  FIX  is  from  0.1  up  to  about  three 
or  four. 

4.  EIGMAX 

Once  the  initial  guess  L  is  calculated,  the  user  has  the  option 
through  IER  of  verifying  that  the  associated  control  matrix  is  stable. 
If  the  user  asks  for  the  stability  check,  the  eigenvalues  of  the  control 
matrix  are  determined.  Each  must  have  its  real  part  more  negative  than 
minus  the  absolute  value  of  EIGMAX.  From  the  discussion  on  cost  matrices, 
theoretically  EIGMAX  could  be  set  to  zero.  Numerically  though,  an 
eigenvalue  that  is  yery   close  to  zero  can  induce  numerical  instability 
in  the  iteration  scheme.  From  experience,  the  suggested  value  is  0.001. 
The  user  can  also  use  FIX  to  modify  the  control  matrix  of  the  initial 
guess  in  an  attempt  to  make  the  real  parts  of  the  eigenvalues  negative 
enough  to  pass  the  stability  check. 

5.  IER 

The  most  informative  parameter  of  the  group  is  IER.  On  return 
from  RICATS,  IER  can  inform  the  user  of  the  validity  of  the  solution. 
When  speicfying  IER  as  in  input  parameter,  the  user  can,  by  bypassing 
various  combinations  of  the  three  checks  available  in  RICATS,  overcome 
some  system  deficiencies,  save  execution  time,  or  obviate  decks  for 
the  subroutines  GMRANK,  ATEIG,  HSBG  (see  note  6  in  comment  cards  at  the 
beginning  of  subroutine  RICATS).  By  bypassing  any  or  all  of  the  three 
checks  (if  the  user  is  fairly  certain  that  the  circumvented  checks 
would  have  been  passed  )  the  user  can  save  execution  time,  although 
the  savings  are  neglibile  except  for  high  order  systems. 
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Care  must  be  exercised  when  not  checking  controllability  and 
or  stability.  For  systems  that  are  uncontrollable,  neglecting  the 
controllability  check  may  lead  to  a  convergent  solution,  if  the  stability 
check  were  passed.  Such  a  solution  could  be  positive  definite  and  yield 
a  stable  final-solution  control  matrix,  but  the  user  should  keep  in 
mind  that  there  may  be  modes  of  the  system  that  cannot  be  controlled. 
For  the  case  of  unstable  control  matrices  corresponding  to  the  initial 
guess,  the  stability  check  should  be  bypassed  only  after  FIX  has  been 
varied  through  the  suggested  range  of  values  without  success.  The 
danger  in  attempting  to  generate  a  solution  for  a  system  that  would  fail 
the  stability  check  is  that  successive  iterations  are  no  longer  bounded 
by  a  lower  positive  [semi]  definite  iteration.  The  iterations  will 
probable  not  converge  to  a  satisfactory  answer.  This  is  the  most 
dangerous  of  the  three  checks  to  bypass,  by  far. 

The  positive  definite  check  of  the  solution  is  included  as  a 
convenience  for  the  user,  to  verify  that  the  final  solution  is  indeed 
positive  definite.  A  note  of  caution  should  be  sounded  when  solving 
high  order  systems.  The  positive  definite  check  is  subject  to  numerical 
problems  when  evaluating  the  high  order  solutions,  and  thus  a  solution 
will  be  flagged  as  not  being  positive  definite,  when  for  all  practical 
purposes  it  is.  Refinements  in  the  solution,  from  introducing  iterative 
routines  for  solving  equation  (16)  (see  subroutine  SIMQ  below),  have 
shown  that  the  difference  between  passing  and  failing  the  check  can 
depend  solely  on  numerically  insignificant  digits  in  a  few  elements  of 
the  solution.  Therefore,  to  give  confidence  to  a  solution  that  has  been 
flagged  as  not  being  positive  definite,  the  user  can:   (1)  determine 
the  eigenvalues  of  the  final  solution,  closed-loop  system  control 
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matrix  and  check  for  all  negative  real  parts;  or  (2)  calculate  the  value 
of  the  optimal  quadratic  cost  function  for  an  arbitrary  set  of  initial 
conditions  using  equation  (15),  and  check  for  a  positive,  finite  result. 
The  positive  definite  check  will  yield  the  same  result  as  the  IBM/SSP 
subroutine  MFSD. 

C.  HIGH  ORDER  SYSTEMS 

As  the  order  of  the  system  increases,  the  problems  due  to  finite 
numerical  calculations  increase  considerably.  One  means  of  trying  to 
maintain  sufficient  accuracy  is  to  have  a  double-precision  version  of 
subroutine  RICATS.  Since  this  would  nearly  double  the  storage  require- 
ments (prohibitive  for  systems  of  order  higher  than  twenty)  it  was  felt 
that  this  was  not  a  feasible  means  of  achieving  the  goal.  The  suggested 
procedure  for  systems  of  higher  order  is  for  the  user  to  introduce  his 
own  dummy  subroutines  MINV  and  SIMQ.  Care  must  be  exercised  when  writing 
your  own  subroutines.  The  variables  passed  to  and  returned  from  your 
subroutine  must  correspond  exactly  to  the  variables  as  handled  by  the 
subroutine  you  are  replacing. 

1 .  Subroutine  MINV 

Using  a  common  statement  from  the  user's  main  program  to  provide 
the  additional  storage  required,  the  user  can  write  a  routine  to  convert 
the  matrix  to  be  inverted  to  double-precision  storage,  invert  the  matrix 
in  double-precision,  then  convert  the  inverted  matrix  back  to  the 
single-precision  mode.  When  this  is  done,  the  inverted  matrix  is 
passed  back  to  the  subroutine  RICATS  as  if  the  IBM/SSP  subroutine  MINV 
had  done  the  inverting.  A  logical  way  to  invert  the  matrix  in  the  double- 
precision  mode  is  to  use  the  IBM/SSP  subroutine  DMINV.  It  might  appear 
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that  this  type  of  routine  would  make  no  noticeable  change  in  the 
execution  of  subroutine  RICATS.  However,  for  some  systems  tested, 
particularly  those  with  singular  initial  guesses  (FIX  =  1.1),  the  number 
of  iterations  required  for  solution  was  halved. 
2.  Subroutine  SIMQ 

Again  using  a  common  statement,  recognizing  that  work  areas  can 
be  declared  single-precision  in  one  subroutine  and  double-precision  in 
another,  the  user  can  write  his  own  simultaneous  linear-equation-solving 
routine.  The  IBM/SSP  subroutine  SIMQ  has  been  found  to  yield  satis- 
factory results  for  a  maximum  of  about  forty-five  equations  (n  =  9). 
For  systems  of  higher  order,  the  solutions  did  converge,  but  they  were 
usually  accompanied  by  an  error  flag  indicating  that  they  were  not 
positive  definite  solutions.  However,  iterative  refinement  of  the 
solution  to  equation  (16)  at  each  iteration  can  yield  error-free 
results.  The  user  can  either  write  his  own  routine  for  the  iterative 
refinement,  or  pass  the  set  of  equations  to  the  IBM/SSP  subroutines 
FACTR  and  RSLMC.  It  is  suggested  that  the  common  statement  also  contain 
a  relative  tolerance  parameter,  not  necessarily  the  same  as  TOLER 
inputed  to  RICATS,  in  either  name  or  magnitude.  As  was  previously 
mentioned,  this  subroutine  technique  was  used  to  remove  error  flags 
from  the  solutions  returned  by  subroutine  RICATS. 
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c 

c      

c 

C         SUBROUTINE  RICATS 
C 

C         P  URPCSE 

C  TO  1TERATIVELY  DETERMINE  THE  UNIQUE,  SYMMETRIC, 

C  POSITIVE-DEFINITE  SOLUTION  P  TO  THE  NONLINEAR, 

C  STEADY-STATE  MATRIX  RICCATI  EQUATION 

C  .                                    -1 

C  P  =  0  =  -  P*A  -  A»*P  -  C«*C  +  P*B*R   *B«*P 

C  OF  OPTIMAL  CONTROL  THEORY.  SEE  NOTE  10. 

C 

C         USAGE 

C  CALL  RICATS  ( A ,6 ,Q ,R ,P ,D , E ,F , V ,L ,L I , MI , I  A, JB , KQ, 

C  1               NTRY,TOLER,FIX, E IGMAX, I ER ,NN ) 

C 

C         DESCRIPTION  OF  PARAMETERS 

C  A       -  GENERAL  IA  BY  IA  INPUT  MATRIX.  SEE  NOTE 

C  3  • 

C  B         GENERAL  IA  BY  JB  INPUT  MATRIX. 

C  Q       -  LOWER  TRIANGULAR  (OR  UPPER  TRIANGULAR, 

C  SEE  KQ)  PART  OF  A  SYMMETRIC,  PCSITVE 

C  SEMIDEFINITE  IA  BY  IA  INPUT  MATRIX. 

C  Q  =  C • *C • 

C  R       -  GENERAL  PART  OF  A  SYMMETRIC,  POSITIVE 

C  DEFINITE  JB  BY  JB  INPUT  MATRIX. 

C  P       -  IA  BY  IA  WORK  MATRIX.  ON  RETURN  P  CON- 

C  TAINS  THE  GENERAL  FORM  OF  THE  SOLUTION. 

C  D       -  MM  BY  MM  WORK  MATRIX.  SEE  NOTE  1. 

C  E       -  MM  BY   1  WORK  VECTOR. 

C  F       -  MZ  BY   1  WORK  VECTOR.  SEE  NOTE  1. 

C  V         MM  BY   1  WORK  VECTOR. 

C  L       -  IA  BY  IA  INTEGER  WORK  MATRIX. 

C  LI      -  IA  BY   1  INTEGER  WORK  VECTOR. 

C  MI      -  IA  BY   1  INTEGER  WORK  VECTOR. 

C  IA      -  ROW  AND  COLUMN  DIMENSION  OF  MATRICES  A, 

C  Q  AND  RETURNED  SOLUTION  P.  ROW 

C  DIMENSION  OF  MATRIX  B. 

C  JB      -  COLUMN  DIMENSION  OF  MATRIX  3.  ROW  AND 

C  COLUMN  DIMENSION  OF  MATRIX  R.  JB  <=  I  A. 

C  KQ      -  INTEGER  INPUT  CONSTANT. 

C  =  0  -  Q  MATRIX  IS  STORED  COLUMNWISE  IN 

C  LOWER  TRIANGULAR  FORM. 

C  =  1  -  Q  MATRIX  IS  STORED  COLUMNWISE  IN 

C  UPPER  TRIANGULAR  FORM  (SEE 

C  I3M/SSP  SUBROUTINE  MSTR). 

C  NTRY    -  MAXIMUM  NO.  OF  ITERATIONS  TO  BE  TRIED. 

C  2  >  NTRY  >  151. 

C  TOLER   -  INPUT  CONSTANT  WHICH  IS  USED  AS  A 

C  RELATIVE  TOLERANCE  FOR  TEST  OF 

C  CONVERGENCE. 

C  FIX     -  INPUT  CONSTANT.  THEORETICALLY  SHOULD  BE 

C  >=  1.0.  SUGGESTED  VALUE  IS  1.1.  SEE 

C  NOTE  7. 

C  EIGMAX  -  MINIMUM  ACCEPTABLE  MAGNITUDE  OF  THE 

C  NEGATIVE  REAL  PARTS  OF  THE  EIGENVALUES 

C  OF  THE  CONTROL  MATRIX  (  A  -  B*V(0)  )  OF 

C  THE  INITIAL  GUESS  V(0).  SUGGESTED  VALUE 

C  OF  EIGMAX  IS  OF  ORDER  0.001  .  SEE  NOTE 

C  5.  EIGMAX  USED  ONLY  FOR  I ER  =  0,1,4,5. 

C  IER     -  INTEGER  INPUT/OUTPUT  PARAMETER. 

C  THE  FOLLOWING  NOTATION  IS  USED  BELOW  : 

C  C.A,B.  -  FOR  NUMERICAL  CONTROLLABILITY 

C  OF  THE  MATRIX  PAIR  (A,B1. 

C  SEE  NOTES  4,6. 

C  S.C.M.  -  FOR  A  STABLE  CONTROL  MATRIX  FOR 

C  THE  INITIAL  GUESS.  SEE  NOTES  5,6. 

C  P.D.S.  -  FOR  A  POSITIVE-DEFINITE  SOLUTION. 

C  SEE  NOTE  8. 

C 
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INPUT 


0 

1 

2 
3 
4 
5 
6 
7 


-  THE  VALUE  OF  IER  DETERMINES  WHICH 
OF  THE  3  ABOVE  PROPERTIES  RICATS 
IS  TO  CHECK  DURING  EXECUTION.  SEE 
NOTE  9. 

S.C.M.  ;  P.D.S. 
;  S.C.M. 


CHECK 
CHECK 
CHECK 
CHECK 
CHECK 
CHECK 
CHECK 
MAKE 


.ATB, 
.A,B 
.  A,B, 
.A,B 


CM. 
CM. 


P.D.S 
P.D.S 


P.D.S. 
NONE  OF  THE  ABOVE  3  CHECKS. 


NN 


OUTPUT  -  ON  R 
FLAG 
ALL  CHE 
SOLUTIO 
ITERATI 
ALL  CHE 
SOLUTIO 
NTRY  IT 
C.A  ,B. 
S.C.M. 
P.D.S. 
CONVERG 
P.D.S. 
NOT  CON 
IMPROPE 
RICATS. 
MORE  TH 
RETURNE 
CAT ING 
ROUTINE 
NOTF  7. 
INTEGER  OUT 
FOR  RICATS 
IS  THE 
BY  SUBR 
RICATS 
IS  THE 
CONTROL 


ETURN  IER  SERVES  AS  AN  ERROR 


=  0 


=  1  - 


=  2  - 

=  3  - 

=  4  - 

=  5  - 

=  6  - 

=  7  - 


FOR 


CKS  REQU 
N  CONVER 
ONS. 

CKS  REQU 
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SUBROUTINE  RICATS  ( A ,B , Q ,R , P ,D , E , F , V,L , LL , M I , I  A, JB , KQ, 
1  NTRY,TOLER,FIXf EIGMAX,IER,NN) 

DIMENSION  A(l) ffBCl)fR(l) «Q(1)«V (lit E(l)yP(l).F(l)9 

1  D  (  1 )  ,  L  {  i  )  ,  L  L  ( 1 )  ,  M  I  (  1  ) 

IF  (    JB.LT.l  .OR.    JB.GT.IA   )  CO  TO  6 

IF  (   IER.LT.O  .OR.   IER.GT.7    )  GO  TO  6 

IF  (  NTRY.LT.2  .OR.  NTRY.GT.150  )  GO  TO  6 

GO  TO  15 
6  IER  =  6 

NN   =  0 

RETURN 
15  CONTINUE 

DIVCK   =  0.1E-30 

NN    =0 

12    =  (IA+D/2 

IAI   =  IA  +  1 

IA2  =  IA  +  2 

MM    =  I  A*  (IA+D/2 

IAS   =  IA*IA 

MMP  =  MM* MM 

EIGMAX  =  -  ABS(EIGMAX) 

IF  (  IA.EQ.l  )  GO  TO  805 

IF  (  IER.GT.3  )  GO  TO  375 
C 

C         IF  CONTROLLABILITY  CHECK  REQUESTED,  FORM  AUGMENTED 
C     CONTROLLABILITY  MATRIX  AND  DETERMINE  ITS  RANK. 

KL  =  IA*JB 

DO  300  1=1, KL 

D( I)  =  B(  I  J 
300    P( I)  =  B(I  ) 

IR  =  KL 

DO  320  M=2,IA 

LI  =  0 

LJ  =  -  IA 

DO  310  J=1,JB 

LJ  =  LJ  +  IA 

DO  310  I  =  1,1  A 

KI  =  I  -  IA 

KJ  =  LJ 

LI  =  LI  +  1 

Ft  LI  )  =  O.OEO 

DO  310  K=1,IA 

KI  =  KI  +  IA 

KJ  =  KJ  +  1 
310    F(LI)  =  F(LI)  +  A(KI)*P(KJ) 

DO  320  1=1, KL 

IR  =  IR  +  1 

D( IR)  =  F(  IJ 
320    PI  I )  =  F{  I  ) 

CALL  GMRANK  ( D , I  A, KL , TOL ER  ,  K) 

IF  {  K.E3. IA  )  GO  TO  375 

NN   =  K 
2  IER  =  2 

RETURN 
C  END  CHECK  ROUTINE 

C 

375  IF  (  KQ.EQ.O  )  GO  TO  35 

DO  2  3  1=1, MM 
20     D(I )  =  Q( I ) 

IR  =  0 

KJ  =  1 

DO  30  J=l, IA 

KJ  =  KJ  +  J  -  1 

KI  =  KJ 

DO  30  I=J,IA 

KI  =  KI  +  I  -  1 

IR  =  IR  +  1 
30     Q( IR)  =  DIKI ) 
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L(KI )  =  IR 
LUR  +  KJ)  =  IR 
IF  (  JB.GT.l  )  GO  TO  45 
D(l)  =  l.OEO/RU) 
GO  TO  55 
KL  =  JB*JB 
DO  50  1=1, KL 
D(I)  =  R(I  ) 

CALL   MINV  (D, JS,DF,LL,MI) 
IR  =  0 
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SAVE  =  FIX*SQRT(  IA)*(MAX  OVER  I  OF  SUM  OVER  J  ABS(A(I,J))) 
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c 
c 


500 
510 

520 


KI  =  -  IA 
DO  100  1=1, IA 
KI  s  KI  +  IAI 
P(KI)  =  P(KI) 


+  SAVE 

SOLVE 

<A+SAVE*I)*S  +  S*(A+SAVE*I )■  =  2*8*6* 
INITIAL  GUESS  V(OJ  =  B« *S {  INVERSE  ) 
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DO  1 
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D(IR 
DO  1 
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D(  IR 
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DO  1 
LJ  = 
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KI  = 
KJ  = 
IR 
P( 
DO 
KI 
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P(  IR 
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IR 
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MIN 

IAS 

-  IA 
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0 
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LJ 
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20  K  = 

KI  + 
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)  =  D 
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IAS 
30  J  = 

LJ  + 
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LJ 

IR  + 
)  =  A 
30  K= 

KI  + 

KJ  + 
)  -  P 

IER. 


PS  (P,D,L,F,IA,MM, IS,MMP, IA2) 

E.O  )  NN  =  1 

l.IAS 

L(I)  ) 

V  (P, I A,DF,LL,MI) 


lflA 

I A 


1,JB 

1 
.GEO 
1,IA 

1 

1 
( IR) 

-  JB 
1,IA 

JB 
ltlA 
IA 

1 
(IRJ 
ltJB 

IA 

1 
(  IR) 
GT.l 


+  B(KI)*P(KJ)*1.0D0 


B(KI)*DIKJ)*1.0D0 

.ND.  IER.LT.4  .OR.  IER.GT.5  )  GO  TO  525 


GUESS  STABILITY  CHECK  REQUESTED,  FORM 
AND  CHECK  ITS  EIGENVALUES. 


C 

c 

525 
140 


IF  INITIAL 
CONTROL  MATRIX 

KI  =  -  IA 

V(2)  =  O.OEO 

DO  500  1=1  ,IA 

KI  =  KI  +  IAI 

V(2)  =  V(2)  +  P(K1  ) 

IF  (  V(2).GT.IA*EIGMAX  )  GO  TO  3 

DO  510  1=1, IAS 

D( I  J  =  P(  I  ) 

CALL   HSBG   (IA,D,IA) 

CALL   ATE1G  ( I  A ,D , V , F , LL , I  A) 

DO  520  1  =  1,  IA 

IF  (  V( I) .GT.EIGMAX  )  GO  TO  3 

CONTINUE 

GO  TO  525 
l  IER  =  3 

KL   =  3 

GO  TO  275 

END  CHECK  ROUTINE 


CONTINUE 

DO  140  1=1, MM 

Q( I)  =  -Q(I) 
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150 
C 

c 
c 


160 


210 


220 


IR  =  0 

LJ  =  IAS  -  JB 
DO  150  J=ltIA 
LJ  =  LJ  +  J8 
DO  150  I  =  J  ,  I  A 
LI  =  ( I-1)*JB  +  IAS 
KJ  =  LJ 
IR  =  IR  +  1 
F( IR)  =  Q(IR) 
KL  =  0 

DO  150  M=1,JB 
KI  =  LI 
KJ  =  KJ  +  1 
DO  150  K=1,JB 
KI  =  KI  +  1 
KL  =  KL  +  1 

FUR)  =  FUR)  -  D(KI)*R(KL)*D(KJ)*1.0D0 

SOLVE 
{A-B*V10))'*V(1)  +  VU)*(A-B*V(0)  )  =  -  Q  -  V(0)*R*V(0) 

CALL  MLIAPS  (  P  ,  D,  L  ,  F, I A,MM , I S, MMP, I A2 ) 

IF  (  IS.NE.O  )  NN  =  NN  +  1 

DO  160  1=1, MM 

IF  (  FU).EQ. O.OEO  )  F(I)  =  DIVCK 

IF  (   ABS(FU)  ).LT.  DIVCK  )  FU)  =   S IGN ( D  I  VCK , F(  I )  ) 

V( I)  =  F(I  ) 

DO  255  M=2,NTRY 

KL  =  0 

IR  =  0 

LJ  =  -  IA 

DO  210  J=1,IA 

LJ  =  LJ  +  IA 

KI  =  0 

DO  210  I  =1,1  A 

KJ  =  LJ 

IR  =  IR  +  1 

PUR)  =  O.OEO 

DO  210  K=i,IA 

KI  =  KI  +  1 

KJ  =  KJ  +  1 

PUR)  =  PUR) 

IR  =  0 

LJ  =  -  IA 

DO  220  J=1,IA 

LJ  =  LJ  +  IA 

DO  220  I=J  ,1  A 

KI  =  ( I-1)*IA 

KJ  =  LJ 

IR  =  IR  +  1 

FUR)  =  Q(IR) 

DO  220  K=1,IA 

KJ  =  KJ  +  1 

KI  =  KI  +  1 

F( IR)  =  F(  IR) 

DO  230  1=1, IAS 

P(I)  =  All)  +  P( I ) 

SOLVE 
VCM+U*  (A+E*V(M))  =  - 


+  E( L(KI ) )*V( L(KJ) 1*1. ODO 


+  V( L(KI ) )*P(KJ)*1.0D0 


V(M+1) 


Q  +  V(M)*E*V(M) 


230 
C 

C  (A+E*V(M)  )  ' 
C 

CALL  MLIAPS  (P , D ,L , F , I  A, MM , IS , MMP,  I A2) 

IF  (  IS.EQ.O  )  GO  TO  235 

NN  =  NN  +  1 

IF     (    NN.GT.IZ    )    GO    TO    7 
235    DO    240    1=1, MM 

IF    {     FU)  .EQ. O.OEO     )     FU)     =    DIVCK 

IF     (        ABS(F(  I)  )  .LT.DIVCK     )     FU)     =       S  IGN  (  D  I  VCK  ,  F  U  )  ) 

IF    (        ABSd.OEO   -    FU)/V(I  )  )  .GT.TOLER    )     GO    TO    245 
240  VII)     =    FU  ) 

GO    TO    265 


33 


245 


250 
255 

265 

270 
275 

280 


DO  25 
IF  ( 
IF  ( 
VU) 
CONTI 
KL  = 
CONTI 
DO  27 
QU) 
IF  ( 
28 
I) 


29 


290 
295 

999 

C 
C 

C     I 


605 

610 
615 

620 


DO 
D( 
IR 

KJ 
DO 
KJ  = 

DO  29 
KI  = 
IR  = 
QUI  ) 
CONTI 
DO  99 
P(I) 
IF  ( 


J  =  I ,MM 
(J).EQ.O.OEO  )  F(J)  =  DIVCK 
ABS(F( J) ) .LT.DIVCK  )  F(J)  = 

F(J) 

UE 


SIGN(DIVCK,F( J) ) 


UE 

1=1, MM 

-Q(  I) 
O.EQ.O  ) 

1=1, MM 

Q(I  ) 


GO  TO  295 


J  = 
J  + 
J 

1  = 
I  + 
R  + 
=  D 
UE 

1  = 

V< 
ER- 


ltlA 

J  -  1 

J,IA 

I  -  1 

1 
(  IR) 

l.IAS 
LCI)  ) 
2*( IER/2) .NE.O 


)  GO  TO  645 


I 

:ACTO 


63  0 

64  0 


C 

C 


IF 
IF 
IF 
GO 
IR 
SAVE 
DO  6 
IR  = 
D(  IR 
IR  = 
I  = 
IR  = 
SAVE 
LI  = 
DO  6 
LI  = 
SAVE 
IF  ( 
IF  ( 
SAVE 
LI  = 
LJ  = 
IZ  = 
M  = 
DO  6 
KI  = 
LJ  = 
KJ  = 
IZ  = 
D(  IZ 
DO  6 
KI  = 
KJ 
D( 
Dl 
I  = 
GO  T 
IER 
RETU 


( 
( 
( 

Tl 
=  1 


POSITIVE-DEFINITE  CHECK  REQUESTED,  ATTEMPT  TO 
THE  SOLUTION  BY  CHOLESKY  SQUARE  ROOT  METHOD. 

LT.DIVCK  )  GO  TO  4 

.2  )  GO  TO  605 

P(3)*P(3)/P(1) .LT.DIVCK  )  GO  TO  4 


(1) 

A.GT 

14)- 

645 


10 
I 

) 
1 

2 
I 

I 

20 
L 


SQRT(Pd)  ) 
J=2,IA 

R  +  IA 

=  P(IR)/SAVE 


IZ 
IZ 


R  +  IAI 

P(  IR) 
R  -  I 

K=2,I 
I  +  1 

=  SAV 

SAVE. 

I  .GE. 

=   SQ 

IR  - 

LI 

IR 

I  +  1 
40  J=M 

LI 

LJ  + 

LJ 

IZ  +  IA 
)  =  P(IZ) 
30  K=2,I 

KI  +  1 

KJ  +  1 
)  =  D(IZ) 
)  =  D(IZ)/SAVE 
I  +  1 
0  615 
=  4  + 
RN 


E  -  DILI )*D(LI ) 
LT.DIVCK  )  GO  TO 
IA  )  GO  TO  645 
RT(SAVE) 

I 


,IA 
IA 


-  D(KI)*D(KJ)*1.0D0 


KL 


END  CHECK  ROUT  INE 
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c 
c 


645  IER  =  KL 
RETURN 
7  IER  =  7 
KL  =  7 
GO  TO  265 

FOR  IA  =  1,  SOLVE  THE  RESULTING  QUADRATIC  EQUATION. 
805  KL  =  0 

IF  (  IER.LT.4  .AND.   ABS (B { 1 ) ) . LT. DI VCK  )  GO  TO  2 

D(l)  =  R(1)/B(1)*B(1) 

D(2)  =  A(1)*D( 1) 

P(l)  =  D(2)  +  SQRT{  D(2)*D(2)  +  Q(l)/D(l)  ) 

IF  (  IER-2*( IER/2) .EQ.O  .AND.  P ( 1) .LT.O.OEO  )  GO  TO 

IER  =  KL 

RETURN 

END 


10 


15 


SUBR 
DIME 
DO  5 
Oil  ) 
IR 


OUTINE  ML  I  APS  IP ,D,L t F , IA, MM, I S , MMP, I A2) 
NSION  P(1),D(1),L(1),F(1) 

1=1, MMP 

=  O.OEO 


=  0 


DO 
LI 
DO 
KJ 
KI 
IR 
DO 

KI  = 
KJ  = 
LJ  = 
D(LJ 
KI  = 
DO  1 
KI  = 
KJ  = 
DO  1 
KJ  = 
D(KJ 
CALL 
RETU 
END 


10 


=  I 

10 
=  J 
=  L 
I 

K 

K 

L 


10 

~) 
5 


K 

K 
5 

K 
) 

S 
RN 


1=1, 

-  I 
J=l, 

-  I 
I 

R  + 
K=l, 
I  + 
J  + 
(KI  ) 
=  D( 

IA 
1  =  1, 
I  + 
I  - 
J  =  l, 
J  + 
=  2. 
IMQ 


IA 
A 
IA 
A 

1 

IA 
IA 
IA 
+ 
LJ) 


(L(KJ }-l)*MM 
+  P( IR) 


IA 

IA2  -  I 

MM 

MM 

MM 

OEO-D(KJ) 

(D,F,MM,  IS) 
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c 

c  , 

c 

C  SUBROUTINE  GMRANK 

C 

C  PURPOSE 

C  TO  DETERMINE  THE  NUMERICAL  ROW  (COLUMN)  RANK  OF 

C  A  GENERAL  MATRIX. 

C 

C  USAGE 

C  CALL  GMRANK  ( D, I  A, KL , TOL ER ,K ) 

C 

C  DESCRIPTION  OF  PARAMETERS 

C  D         GENERAL  INPUT  MATRIX  (DESTROYED). 

C  IA      -  ROW  DIMENSION  OF  MATRIX  D. 

C  KL        COLUMN  DIMENSION  OF  MATRIX  D. 

C  TOLER   -  INPUT  CONSTANT  USED  AS  A  RELATIVE 

C  TOLERANCE  FOR  LOSS  OF  SIGNIFICANCE. 

C  K       -  ON  RETURN,  K  =  RANK  OF  MATRIX  D. 

C 

C  REMARKS 

C  NONE 

C 

C  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 

C  NONE 

C 

C  REFERENCES 

C  PENNINGTON, RALPH  H.,  "INTRODUCTORY  COMPUTER 

C  METHODS  AND  NUMERICAL  ANALYSIS",  PUBLISHED  BY 

C  MACMILLAN  COMPANY,  NEW  YORK,  1965. 

C 

C      

c 


SUBROUTINE  GMRANK  ( D , I  A , KL , TOL ER ,K ) 

DIMENSION  Dll) 

IAI  =  I A  +  1 

NN  =  MINOl  IA,KL) 

IR  =  -  IA 

K  =  1 
325  CONTINUE 

IR  =  I R  +  IAI 

SAVE  =   ABS(DUR)) 

LI  =  K 

LJ  =  K 

KJ  =  K  -  1 

KI  =  IR  -  K 

DO  330  J=K,KL 

KI  =  KI  +  KJ 

DO  330  I=K,IA 

KI  =  KI  +  1 

IF  (   ABS(D(KI) ) .LE.SAVE  )  GO  TO  330 

LI  =  I 

LJ  =  J 

SAVE  =   ABSID(KI)) 
330    CONTINUE 

IF  (  SAVE.LT.0.1E-35  )  GO  TO  365 

IF  (  K.GE.NN  )  GO  TO  375 

IF  (  LI.LE.K  )  GO  TO  345 

KI  =  IK  -  IA 

KJ  =  KI  +  LI  -  K 

DO  340  J=K,KL 

KJ  -    KJ  +  IA 

KI  =  KI  +  IA 

SAVE   =  D(KI ) 

D(KI )  =  D(KJ) 
340    D(KJ)  =  SAVE 
345  CONTINUE 

IF  (  LJ.LE.K  )  GO  TO  355 
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350 
355 


360 


365 
375 


KI 

KJ 

DO 

KI 

KJ  = 

SAVE 

D(KI 

D(KJ 


=  IR  -  1 

=  KI  +  (LJ-K)*IA 

350  I=K,IA 

=  KI  +  1 

=  KJ  +  1 
=  D(KI) 
)  =  DtKJJ 
)  =  SAVE 


CONTINUE 

LI  =  K  +  1 

LJ  =  IR 

KJ  =  IR  +  I A  -  K 

DO  360  J=LI ,KL 

LJ  =  LJ  +  IA 

KI  =  IR 

KJ  =  KJ  +  K 

D(LJ )     =  D(LJ)/D( IR) 

DO  360  I=LI, IA 

KI  =  KI  +  1 

KJ  =  KJ  +  1 

SAVE   =  D(LJ)*D(KI ) 

D(KJ)  =  D(KJ)  -  SAVE 

IF  (   ABS(D(KJ ) ) .GE.TOLER* 

D(KJ)  =  O.OEO 

CONTINUE 

K  =  K  +  1 

GO  TO  325 

CONTINUE 

K  =  K  —  1 

CONTINUE 

RETURN 

END 


ABS(SAVE)  )  GO  TO  360 
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